#ifndef FUNCTION_H_
#define FUNCTION_H_

#include "VarArray.hpp"
#include "defs.hpp"
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace Eigen;
using namespace std;
void GridGenerate(VarArray<Point, 2, 0> &grid, int np_in_xdir, int np_in_ydir);

void OutPutGridInPlot3d(const VarArray<Point, 2, 0> &grid);

void InitializeFlow(VarArray<Variable, 2, 0> &flow);

void ThomasAlgorithm(VarArray<double, 2, 0> &coefficient_matrix, double *value,
                     int begin, int end);

void InitializeNextLine(VarArray<Variable, 2, 0> &flow, int column);

struct Feta {
  double s_{6.0};
  double sx_{4.0};
  double y_max_ = 1.0;
  int np_in_ydir_, np_in_xdir_;
  double delty, deltx;
  Feta(double s, double ymax, int np_in_ydir, double sx, int np_in_xdir)
      : s_{s}, y_max_{ymax}, np_in_ydir_{np_in_ydir}, sx_{sx}, np_in_xdir_{
                                                                   np_in_xdir} {
    delty = 1.0 / (np_in_ydir_ - 1);
    deltx = 1.0 / (np_in_xdir_ - 1);
  }
  double operator()(int j) { return y_max_ * sinh(j * delty * s_) / sinh(s_); }
  double f1(double j) { return y_max_ / sinh(s_) * s_ * cosh(s_ * j * delty); }
  double f2(double j) {
    return y_max_ / sinh(s_) * s_ * s_ * sinh(s_ * j * delty);
  }

  double getX(int i) { return 2.0 * sinh(i * deltx * sx_) / sinh(sx_); }
};
/**
class CalcUVsecondOrder {
private:
  VarArray<Variable, 2, 0> &flow;
  const VarArray<Point, 2, 0> &grid;
  int np_in_xdir, np_in_ydir;
  Feta feta;
  VarArray<double, 2, 0> coef;
  double re_num;
  const double delt_eta;
  vector<double> fx_u;
  vector<double> nut;

public:
  CalcUVsecondOrder(VarArray<Variable, 2, 0> &flow_,
                    const VarArray<Point, 2, 0> &grid_, double re_num_)
      : flow{flow_}, grid{grid_}, np_in_xdir{flow.Dim(0)}, np_in_ydir{flow.Dim(
                                                               1)},
        feta(10, 1.0, np_in_ydir, 7, np_in_xdir), delt_eta{1.0 /
                                                           (np_in_ydir - 1)},
        fx_u(np_in_ydir), nut(np_in_ydir) {
    coef.ReConstruct(3, np_in_ydir);
    re_num = re_num_;
  }

  void CalcNut(int column) {
    static constexpr double ka{0.4};
    nut[0] = 0.0;
    nut[np_in_ydir - 1] = 0.0;
    for (int i = 1; i < np_in_ydir - 1; ++i) {
      nut[i] = (ka * grid(column, i).y) * (ka * grid(column, i).y) *
               fabs((flow(column, i + 1).u - flow(column, i - 1).u) * 0.5 /
                    delt_eta / feta.f1(i));
      nut[i] *= re_num;
    }
  }

  void CalcCoef1(int column) {

    double rdelt_eta = 1.0 / delt_eta;
    double C = 1.0 / re_num;

    for (int j = 1; j < np_in_ydir - 1; ++j) {
      double rf1 = 1.0 / feta.f1(j);
      double A = flow(column, j).v * rf1 +
                 C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
                 C * (nut[j + 1] - nut[j - 1]) * 0.5 * rdelt_eta * rf1 * rf1;
      A *= 0.5;
      double B = -C * rf1 * rf1 * (1 + nut[j]);
      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

      coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      coef(1, j) = flow(column, j).u * rdeltx - 2.0 * B * rdelt_eta * rdelt_eta;
      coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      //   cout << coef(0, j) << "  " << coef(1, j) << "  " << coef(2, j) <<
      //   endl;
    }

    coef(0, 1) = 0.0;
    // 0 梯度边界条件

    int j = np_in_ydir - 1;
    double rf1 = 1.0 / feta.f1(j);
    double A = flow(column, j).v * rf1 +
               C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
               C * (nut[j] - nut[j - 1]) * rdelt_eta * rf1 * rf1;
    A *= 0.5;
    double B = -C * rf1 * rf1 * (1 + nut[j]);
    double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

    coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
    coef(1, j) = flow(column, j).u * rdeltx - 2.0 * B * rdelt_eta * rdelt_eta;
    coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;

    coef(1, np_in_ydir - 1) += coef(2, np_in_ydir - 1);
    coef(2, np_in_ydir - 1) = 0.0;
  }

  void CalcCoef2(int column) {

    double rdelt_eta = 1.0 / delt_eta;
    double C = 1.0 / re_num;

    for (int j = 1; j < np_in_ydir - 1; ++j) {
      double rf1 = 1.0 / feta.f1(j);
      double A = flow(column, j).v * rf1 +
                 C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
                 C * (nut[j + 1] - nut[j - 1]) * 0.5 * rdelt_eta * rf1 * rf1;
      A *= 0.5;
      double B = -C * rf1 * rf1 * (1 + nut[j]);
      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 2, j).x);

      coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      coef(1, j) =
          flow(column, j).u * rdeltx * 3.0 - 2.0 * B * rdelt_eta * rdelt_eta;
      coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      //   cout << coef(0, j) << "  " << coef(1, j) << "  " << coef(2, j) <<
      //   endl;
    }

    coef(0, 1) = 0.0;
    // 0 梯度边界条件

    int j = np_in_ydir - 1;
    double rf1 = 1.0 / feta.f1(j);
    double A = flow(column, j).v * rf1 +
               C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
               C * (nut[j] - nut[j - 1]) * rdelt_eta * rf1 * rf1;
    A *= 0.5;
    double B = -C * rf1 * rf1 * (1 + nut[j]);
    double rdeltx = 1.0 / (grid(column, j).x - grid(column - 2, j).x);

    coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
    coef(1, j) =
        flow(column, j).u * rdeltx * 3.0 - 2.0 * B * rdelt_eta * rdelt_eta;
    coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;

    coef(1, np_in_ydir - 1) += coef(2, np_in_ydir - 1);
    coef(2, np_in_ydir - 1) = 0.0;
  }

  void CalcFx1(int column) {
    for (int j = 1; j < np_in_ydir; ++j) {

      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);
      fx_u[j] = flow(column, j).u * flow(column - 1, j).u * rdeltx;
      //   cout << fx_u[j] << endl;
    }
  }

  void CalcFx2(int column) {
    for (int j = 1; j < np_in_ydir; ++j) {

      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 2, j).x);
      fx_u[j] = -flow(column, j).u *
                (flow(column - 2, j).u - 4.0 * flow(column - 1, j).u) * rdeltx;
      //   cout << fx_u[j] << endl;
    }
  }

  void IterateForUV(int column) {

    vector<double> old_u(np_in_ydir);
    for (int j = 1; j < np_in_ydir; ++j)
      old_u[j] = flow(column, j).u;

    int count{0};
    double tmpu, tmpv, res_u, res_v, res_u1, res_v1;
    do {
      CalcNut(column);
      if (column == 1) {
        CalcCoef1(column);
        CalcFx1(column);
      } else {
        CalcCoef2(column);
        CalcFx2(column);
      }

      ThomasAlgorithm(coef, fx_u.data(), 1, fx_u.size());
      //   cout << "hello2" << endl;

      for (int j = 1; j < np_in_ydir; ++j) {
        // cout << fx_u[j] << endl;
        flow(column, j).u = fx_u[j];
      }

      res_u = Diffence(old_u, fx_u);
      res_v = UpdateV(column);
      if (count == 1) {
        res_u1 = res_u + 1.0e-30;
        res_v1 = res_v + 1.0e-30;
      }
      tmpu = log10(res_u / res_u1);
      tmpv = log10(res_v / res_v1);
      ++count;
    } while (tmpu > -12.0 && tmpv > -12.0 && res_u > 1.0e-14 &&
             res_v > 1.0e-14 && count < 1000);

    // cout << scientific << setprecision(14);
    // for (int j = 0; j < np_in_ydir; ++j)
    //   cout << flow(1, j).u << endl;
    // exit(EXIT_FAILURE);
    cout << "count = " << count << endl;
  }

  double Diffence(const vector<double> &v1, const vector<double> &v2) {
    double res{0.0};
    for (int j = 1; j < v1.size(); ++j)
      res += (v1[j] - v2[j]) * (v1[j] - v2[j]);
    return sqrt(res);
  }

  double UpdateV(int column) {

    vector<double> old_v(np_in_ydir);
    for (int j = 1; j < np_in_ydir; ++j)
      old_v[j] = flow(column, j).v;

    double res{0.0};
    for (int j = 1; j < np_in_ydir; ++j) {
      double deltx = (grid(column, j).x - grid(column - 1, j).x);
      flow(column, j).v =
          flow(column, j - 1).v -
          delt_eta * feta.f1(j - 0.5) / deltx *
              (0.5 * (flow(column, j).u + flow(column, j - 1).u) -
               0.5 * (flow(column - 1, j).u + flow(column - 1, j - 1).u));
      // flow(column, j).v = flow(column, j - 1).v -
      //                     delt_eta * feta.f1(j) / deltx *
      //                         ((flow(column, j).u - flow(column - 1, j).u));

      //   cout << flow(column, j).v << endl;
      res += (old_v[j] - flow(column, j).v) * (old_v[j] - flow(column, j).v);
    }
    return sqrt(res);
  }
    double UpdateV2(int column) {

      int npoint_in_ydir{flow.Dim(1)};

      vector<double> old_v(npoint_in_ydir);
      for (int j = 1; j < npoint_in_ydir; ++j)
        old_v[j] = flow(column, j).v;

      //   for (int j = 0; j < npoint_in_ydir; ++j)
      //     cout << flow(column, j).v << endl;
      //   cout << "end v" << endl;

      int size = npoint_in_ydir - 1;
      MatrixXd coef = MatrixXd::Zero(size, size);
      coef(0, 0) = -3;
      coef(0, 1) = 6;
      coef(0, 2) = -1;

      for (int j = 1; j < size - 2; ++j) {
        coef(j, j - 1) = -2;
        coef(j, j) = -3;
        coef(j, j + 1) = 6;
        coef(j, j + 2) = -1;
      }

      coef(size - 2, size - 3) = -2;
      coef(size - 2, size - 2) = -3;
      coef(size - 2, size - 1) = 5;

      coef(size - 1, size - 2) = -2;
      coef(size - 1, size - 1) = 2;

      VectorXd fxb(size);
      for (int j = 0; j < size; ++j) {
        int jp1 = j + 1;

        const double deltx = (grid(column, jp1).x - grid(column - 1, jp1).x);
        fxb(j) = 6.0 * delt_eta / deltx *
                 (flow(column, jp1).u - flow(column - 1, jp1).u);
      }
      fxb *= -1.0;

      VectorXd new_v = coef.colPivHouseholderQr().solve(fxb);

      // update
      for (int j = 0; j < size; ++j)
        flow(column, j + 1).v = new_v(j);

      // for (int j = 0; j < npoint_in_ydir; ++j)
      //   cout << flow(column, j).v << endl;
      // cout << "end v" << endl;
      // exit(EXIT_FAILURE);

      double res_v{0.0};
      for (int j = 1; j < npoint_in_ydir; ++j)
        res_v += (flow(column, j).v - old_v[j]) * (flow(column, j).v -
    old_v[j]); return sqrt(res_v);
    }
}; **/

class CalcUVsecondOrder {
private:
  VarArray<Variable, 2, 0> &flow;
  const VarArray<Point, 2, 0> &grid;
  int np_in_xdir, np_in_ydir;
  Feta feta;
  VarArray<double, 2, 0> coef;
  double re_num;
  const double delt_eta;
  vector<double> fx_u;

public:
  vector<double> nut;
  CalcUVsecondOrder(VarArray<Variable, 2, 0> &flow_,
                    const VarArray<Point, 2, 0> &grid_, double re_num_)
      : flow{flow_}, grid{grid_}, np_in_xdir{flow.Dim(0)}, np_in_ydir{flow.Dim(
                                                               1)},
        feta(10, 1.0, np_in_ydir, 7, np_in_xdir), delt_eta{1.0 /
                                                           (np_in_ydir - 1)},
        fx_u(np_in_ydir), nut(np_in_ydir) {
    coef.ReConstruct(3, np_in_ydir);
    re_num = re_num_;
  }

  void CalcNut(int column) {
    static constexpr double ka{0.4};
    nut[0] = 0.0;
    nut[np_in_ydir - 1] = 0.0;
    for (int i = 1; i < np_in_ydir - 1; ++i) {
      nut[i] = (ka * grid(column, i).y) * (ka * grid(column, i).y) *
               fabs((flow(column, i + 1).u - flow(column, i - 1).u) * 0.5 /
                    delt_eta / feta.f1(i));
      nut[i] *= re_num;
    }
  }

  void CalcCoef(int column) {

    double rdelt_eta = 1.0 / delt_eta;
    double C = 1.0 / re_num;

    for (int j = 1; j < np_in_ydir - 1; ++j) {
      double rf1 = 1.0 / feta.f1(j);
      double A = flow(column, j).v * rf1 +
                 C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
                 C * (nut[j + 1] - nut[j - 1]) * 0.5 * rdelt_eta * rf1 * rf1;
      A *= 0.5;
      double B = -C * rf1 * rf1 * (1 + nut[j]);
      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

      coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      coef(1, j) = flow(column, j).u * rdeltx - 2.0 * B * rdelt_eta * rdelt_eta;
      coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;
      //   cout << coef(0, j) << "  " << coef(1, j) << "  " << coef(2, j) <<
      //   endl;
    }

    coef(0, 1) = 0.0;
    // 0 梯度边界条件

    int j = np_in_ydir - 1;
    double rf1 = 1.0 / feta.f1(j);
    double A = flow(column, j).v * rf1 +
               C * feta.f2(j) * rf1 * rf1 * rf1 * (1 + nut[j]) -
               C * (nut[j] - nut[j - 1]) * rdelt_eta * rf1 * rf1;
    A *= 0.5;
    double B = -C * rf1 * rf1 * (1 + nut[j]);
    double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

    coef(0, j) = -A * rdelt_eta + B * rdelt_eta * rdelt_eta;
    coef(1, j) = flow(column, j).u * rdeltx - 2.0 * B * rdelt_eta * rdelt_eta;
    coef(2, j) = A * rdelt_eta + B * rdelt_eta * rdelt_eta;

    coef(1, np_in_ydir - 1) += coef(2, np_in_ydir - 1);
    coef(2, np_in_ydir - 1) = 0.0;
  }

  void CalcFx(int column) {
    for (int j = 1; j < np_in_ydir; ++j) {

      double rdeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);
      fx_u[j] = flow(column, j).u * flow(column - 1, j).u * rdeltx;
      //   cout << fx_u[j] << endl;
    }
  }

  void IterateForUV(int column) {

    vector<double> old_u(np_in_ydir);
    for (int j = 1; j < np_in_ydir; ++j)
      old_u[j] = flow(column, j).u;

    int count{0};
    double tmpu, tmpv, res_u, res_v, res_u1, res_v1;
    do {
      CalcNut(column);
      CalcCoef(column);
      CalcFx(column);

      ThomasAlgorithm(coef, fx_u.data(), 1, fx_u.size());
      //   cout << "hello2" << endl;

      for (int j = 1; j < np_in_ydir; ++j) {
        // cout << fx_u[j] << endl;
        flow(column, j).u = fx_u[j];
      }

      res_u = Diffence(old_u, fx_u);
      res_v = UpdateV(column);
      if (count == 1) {
        res_u1 = res_u + 1.0e-30;
        res_v1 = res_v + 1.0e-30;
      }
      tmpu = log10(res_u / res_u1);
      tmpv = log10(res_v / res_v1);
      ++count;
    } while (tmpu > -12.0 && tmpv > -12.0 && res_u > 1.0e-14 &&
             res_v > 1.0e-14 && count < 300);

    // cout << scientific << setprecision(14);
    // for (int j = 0; j < np_in_ydir; ++j)
    //   cout << flow(1, j).u << endl;
    // exit(EXIT_FAILURE);
    cout << "count = " << count << endl;
  }

  double Diffence(const vector<double> &v1, const vector<double> &v2) {
    double res{0.0};
    for (int j = 1; j < v1.size(); ++j)
      res += (v1[j] - v2[j]) * (v1[j] - v2[j]);
    return sqrt(res);
  }

  double UpdateV(int column) {

    vector<double> old_v(np_in_ydir);
    for (int j = 1; j < np_in_ydir; ++j)
      old_v[j] = flow(column, j).v;

    double res{0.0};
    for (int j = 1; j < np_in_ydir; ++j) {
      double deltx = (grid(column, j).x - grid(column - 1, j).x);
      flow(column, j).v =
          flow(column, j - 1).v -
          delt_eta * feta.f1(j - 0.5) / deltx *
              (0.5 * (flow(column, j).u + flow(column, j - 1).u) -
               0.5 * (flow(column - 1, j).u + flow(column - 1, j - 1).u));
      // flow(column, j).v = flow(column, j - 1).v -
      //                     delt_eta * feta.f1(j) / deltx *
      //                         ((flow(column, j).u - flow(column - 1, j).u));

      //   cout << flow(column, j).v << endl;
      res += (old_v[j] - flow(column, j).v) * (old_v[j] - flow(column, j).v);
    }
    return sqrt(res);
  }
};

#endif